#include "base.h"

int main(int argc,char **argv)
{

//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@///
          //////////////////
          // Input   ///////
          //////////////////

          double beta, Uplaq;
          char gauge_fixing[1024];
          char gauge_action[1024];
	  int NOa2;
	  char input_momdat[1024];
          char momentum_order[1024];


          ////// Parsing command line

          for(int i=1; i<argc; i++) {
            if(strncmp(argv[i],"-gauge_fixing",13)==0)  sscanf(argv[i],"-gauge_fixing=%s",gauge_fixing);
            if(strncmp(argv[i],"-gauge_action",13)==0)  sscanf(argv[i],"-gauge_action=%s",gauge_action);
            if(strncmp(argv[i],"-T",2)==0)  sscanf(argv[i],"-T=%d",&T);
            if(strncmp(argv[i],"-L",2)==0)  sscanf(argv[i],"-L=%d",&L);
            if(strncmp(argv[i],"-beta",5)==0)  sscanf(argv[i],"-beta=%lf",&beta);
            if(strncmp(argv[i],"-Uplaq",6)==0)  sscanf(argv[i],"-Uplaq=%lf",&Uplaq);
            if(strncmp(argv[i],"-NOa2",5)==0)  sscanf(argv[i],"-NOa2=%d",&NOa2);
	    if(strncmp(argv[i],"-momdat",7)==0)  sscanf(argv[i],"-momdat=%s",input_momdat);
            if(strncmp(argv[i],"-momentum_order",15)==0)  sscanf(argv[i],"-momentum_order=%s",momentum_order);


          }
          printf("==============================\n");
          printf("beta=%0.3f Uplaq=%0.5f\n",beta,Uplaq);
          printf("==============================");


  /////////////////////////
  //// Momenta set
  //////////////////////////////
  char momentafile[1024]; sprintf(momentafile,"%s",input_momdat);
  read_momenta_data(momentafile);
  double *nn=(double*)malloc(sizeof(double)*(Nmom));
  compute_nn(momentum_order,nn);
  Momentum(momentum_order);

  char Oa2print[1024];


     //////////////// substracting terms     
     double ***CorOa2_Zq=malloc(sizeof(double**)*(NOa2));
     for (int iOa2=0; iOa2<NOa2; iOa2++) CorOa2_Zq[iOa2]=malloc(sizeof(double*)*(Nmom));
     for (int iOa2=0; iOa2<NOa2; iOa2++)  for (int imom=0; imom<Nmom; imom++) CorOa2_Zq[iOa2][imom]=malloc(sizeof(double)*(1));
     double ***CorOa2_Lambda=malloc(sizeof(double**)*(NOa2));
     for (int iOa2=0; iOa2<NOa2; iOa2++) CorOa2_Lambda[iOa2]=malloc(sizeof(double*)*(Nmom));
     for (int iOa2=0; iOa2<NOa2; iOa2++)  for (int imom=0; imom<Nmom; imom++) CorOa2_Lambda[iOa2][imom]=malloc(sizeof(double)*(5));

     //iOa2=0 -> uncorrected
     for (int imom=0; imom<Nmom; imom++){
     CorOa2_Zq[0][imom][0]=0.;
     for (int igamma=0; igamma<5 ; igamma++) CorOa2_Lambda[0][imom][igamma]=0.;
     }
     //iOa2=1 corrected
     //iOa2=2 corrected boosted
     double g2;
     for (int iOa2=1; iOa2<NOa2; iOa2++)
     {
        for (int imom=0; imom<Nmom; imom++)
        {
          if(iOa2==1) g2=6.0/beta;
          if(iOa2==2) g2=(6.0/beta)/Uplaq;
          correction_Oa2g2(CorOa2_Zq[iOa2][imom],CorOa2_Lambda[iOa2][imom],g2,SinP2[imom],SinP4[imom],nn[imom],gauge_action,gauge_fixing);
        }

     }
    
    printf("====================================================\n\n\n");

        /////////////////
        // Write output
        ////////////////

   char output[1024];


   for (int iOa2=0; iOa2<NOa2; iOa2++)
   {
    whichOa2(iOa2,Oa2print);
         sprintf(output,"subs_%0.0f_L%dT%d_Zq_%s.out",beta*100,L,T,Oa2print);
	 printf("Output in %s\n",output);
	 FILE *fout=NULL;
	 fout=fopen(output,"w");
                if(fout == (FILE *)NULL) {
                        fprintf(stderr,"Unable to open file %s for writing\n",output);
                        return EXIT_FAILURE;
                 }

	 for (int imom=0; imom<Nmom; imom++)
	{
		fprintf(fout,"%8.10f %8.10f \n",SinP2[imom],CorOa2_Zq[iOa2][imom][0]);

	}
	fclose(fout);

   }

   char gammaprint[8];
   for (int iOa2=0; iOa2<NOa2; iOa2++)
   {
    whichOa2(iOa2,Oa2print);
    for (int igamma=0; igamma<5; igamma++)
    {
	whichGammaProm(igamma,gammaprint);
         sprintf(output,"subs_%0.0f_L%dT%d_Lambda%s_%s.out",beta*100,L,T,gammaprint,Oa2print);
	 printf("Output in %s\n",output);
         FILE *fout=NULL;
         fout=fopen(output,"w");
                if(fout == (FILE *)NULL) {
                        fprintf(stderr,"Unable to open file %s for writing\n",output);
                        return EXIT_FAILURE;
                 }

         for (int imom=0; imom<Nmom; imom++)
        {
                fprintf(fout,"%8.10f %8.10f \n",SinP2[imom],CorOa2_Lambda[iOa2][imom][igamma]);

        }
        fclose(fout);

   }
   }

    return 0;

}

